clear all
set more off
cd /Users/zimaoxiao/Dropbox/CC_Yield_Predict/replication_package/  // Navigate to replication folder on your own machine
global weather "gdd0_24 gdd24_26 gdd26 prec prec2"

qui {
* ============================
* 	  Linear trends spec.
* ============================
use if inrange(year,1950,2020) using data/regression/reg_2050_prediction_2spline, clear
	
* apply trends and yield data filter
merge m:1 county_fips using data/yield/filter, keep(3) nogenerate
	
* regression
reghdfe lncornyield $weather state_fips#c.trend [aweight=corn_areaHarv], absorb(county_fips) vce(cluster county_fips StateXYear)

* save the estimation results
tempfile corn_st_lin
parmest, label saving(`"`corn_st_lin'"',replace) format(p %8.2f) stars(0.1 0.05 0.01)

* extract temperature coefficients		
use parm estimate using `corn_st_lin', clear
keep if _n<4
foreach i in gdd0_24 gdd24_26 gdd26 {
	gen coef_`i' = estimate if parm == "`i'"
	sum coef_`i'
	local k = r(mean)
	replace coef_`i' = `k'
}	
keep if _n==1
drop parm estimate
tempfile corntempcoef
save `corntempcoef', replace

* extract county linear trends coefficients		
use parm estimate using `corn_st_lin', clear		
drop if _n<6
gen quad = ustrright(parm,1)
preserve
drop if quad != "d"
drop quad
destring parm, generate(state_fips) ignore("b.state_fips#c.trend")
drop parm			
rename estimate stlin
tempfile cornlintrend
save `cornlintrend', replace
restore
						
* put coefficients and calculated climate change together
use county_fips CC* using data/regression/reg_2050_prediction_2spline, clear
bys county_fips: keep if _n==1
gen state_fips = floor(county_fips/1000)
merge m:1 state_fips using `cornlintrend', keep(3) noreport nogenerate
cross using `corntempcoef'

foreach k in ssp245 ssp585 {
	preserve
	* calculate climate change impacts by county
	bys county_fips: gen ccimpact_corn_`k' = (CC_gdd0_24_`k'*coef_gdd0_24) + (CC_gdd24_26_`k'*coef_gdd24_26) + (CC_gdd26_`k'*coef_gdd26)
					
	* calculate overall impacts by county (climate change + tech effects)
	* 1950 as the 1st year, t=1; 2020 as the 71st year, t=71; 2050 as the 101th year, t=101
	bys county_fips: gen trends_corn_`k' = (CC_gdd0_24_`k'*coef_gdd0_24) + (CC_gdd24_26_`k'*coef_gdd24_26) + (CC_gdd26_`k'*coef_gdd26) + (stlin*30)
	
	keep county_fips ccimpact_corn_`k' trends_corn_`k' 	
	tempfile corn_cccoef_`k'
	save `corn_cccoef_`k'', replace
	restore	
}
	
* merge in regression coefficients, calculate predicted 2020 yields
use data/yield/filter, clear
merge 1:1 county_fips using `corn_cccoef_ssp245', keep(3) noreport nogenerate
merge 1:1 county_fips using `corn_cccoef_ssp585', keep(3) noreport nogenerate
	
* ssp245
* predicted 2050 yields: climate change
bys county_fips: gen cc2050_lin_245 = cornyield_2020*(1+ccimpact_corn_ssp245)

* predicted 2050 yields: climate change + tech effects
bys county_fips: gen trends2050_lin_245 = cornyield_2020*(1+trends_corn_ssp245)

* ssp585
* predicted 2050 yields: climate change
bys county_fips: gen cc2050_lin_585 = cornyield_2020*(1+ccimpact_corn_ssp585)

* predicted 2050 yields: climate change + tech effects
bys county_fips: gen trends2050_lin_585 = cornyield_2020*(1+trends_corn_ssp585)
		
* Table S3: rows (1), (2), (4), (6), (8)
rename (cornyield_2020) (yield2020)
keep county_fips yield2020 cc2050_lin_245 cc2050_lin_585 trends2050_lin_245 trends2050_lin_585
tempfile twopiecewise_2050_lin
save `twopiecewise_2050_lin', replace
		
* ============================
* 	 Quadratic trends spec.
* ============================
use if inrange(year,1950,2020) using data/regression/reg_2050_prediction_2spline, clear
	
* apply trends and yield data filter
merge m:1 county_fips using data/yield/filter, keep(3) nogenerate
				
* regression
reghdfe lncornyield $weather state_fips#c.trend state_fips#c.trend2 [aweight=corn_areaHarv], absorb(county_fips) vce(cluster county_fips StateXYear)

* save the estimation results
tempfile corn_st_quad
parmest, label saving(`"`corn_st_quad'"',replace) format(p %8.2f) stars(0.1 0.05 0.01)
											
* extract temperature coefficients		
use parm estimate using `corn_st_quad', clear
keep if _n<4
foreach i in gdd0_24 gdd24_26 gdd26 {
	gen coef_`i' = estimate if parm == "`i'"
	sum coef_`i'
	local k = r(mean)
	replace coef_`i' = `k'
}	
keep if _n==1
drop parm estimate
tempfile corntempcoef
save `corntempcoef', replace
						
* extract county linear trends coefficients		
use parm estimate using `corn_st_quad', clear
drop if _n<6
gen quad = ustrright(parm,1)
preserve
	* first-order
	drop if quad != "d"
	drop quad
	destring parm, generate(state_fips) ignore("b.state_fips#c.trend")
	drop parm			
	rename estimate stlin
	tempfile cornlintrend
	save `cornlintrend', replace
restore
	* second-order
	keep if quad == "2"
	destring parm, generate(state_fips) ignore("b.state_fips#c.trend")
	gen state_fips_1 = floor(state_fips/10)
	replace state_fips = state_fips_1
	drop state_fips_1
	drop parm quad
	rename estimate stquad
	tempfile cornquadtrend
	save `cornquadtrend', replace
	
* put coefficients and calculated climate change together
use county_fips CC* using data/regression/reg_2050_prediction_2spline, clear
bys county_fips: keep if _n == 1
gen state_fips = floor(county_fips/1000)
merge m:1 state_fips using `cornlintrend', keep(3) noreport nogenerate
merge m:1 state_fips using `cornquadtrend', assert(3) keep(3) noreport nogenerate
cross using `corntempcoef'
	
foreach k in ssp245 ssp585 {
	preserve
	* calculate climate change impacts by county
	bys county_fips: gen ccimpact_corn_`k' = (CC_gdd0_24_`k'*coef_gdd0_24) + (CC_gdd24_26_`k'*coef_gdd24_26) + (CC_gdd26_`k'*coef_gdd26)
		
	* calculate overall impacts by county (climate change + tech effects)
	* 1950 as the 1st year, t=1; 2020 as the 71st year, t=71; 2050 as the 101th year, t=101
	bys county_fips: gen trends_corn_`k' = (CC_gdd0_24_`k'*coef_gdd0_24) + (CC_gdd24_26_`k'*coef_gdd24_26) + (CC_gdd26_`k'*coef_gdd26) + (stlin*30) + (stquad*5160)

	keep county_fips ccimpact_corn_`k' trends_corn_`k' 
	tempfile corn_cccoef_`k'
	save `corn_cccoef_`k'', replace
	restore	
}
	
* merge in regression coefficients, calculate predicted 2050 yields
* calculate predicted 2050 yield
use data/yield/filter, clear
merge 1:1 county_fips using `corn_cccoef_ssp245', keep(3) noreport nogenerate
merge 1:1 county_fips using `corn_cccoef_ssp585', keep(3) noreport nogenerate
	
* ssp245
* predicted 2050 yields: under climate change
bys county_fips: gen cc2050_quad_245 = cornyield_2020*(1+ccimpact_corn_ssp245)

* predicted 2050 yields: under climate change + tech effects
bys county_fips: gen trends2050_quad_245 = cornyield_2020*(1+trends_corn_ssp245)

* ssp585
* predicted 2050 yields: under climate change
bys county_fips: gen cc2050_quad_585 = cornyield_2020*(1+ccimpact_corn_ssp585)

* predicted 2050 yields: under climate change + tech effects
bys county_fips: gen trends2050_quad_585 = cornyield_2020*(1+trends_corn_ssp585)

* Table S3: rows (3), (5), (7), (9)	
rename (cornyield_2020) (yield2020)
keep county_fips cc2050_quad_245 cc2050_quad_585 trends2050_quad_245 trends2050_quad_585
tempfile twopiecewise_2050_quad
save `twopiecewise_2050_quad', replace

* ============================
* 	         Output
* ============================
use `twopiecewise_2050_lin', clear
merge 1:1 county_fips using `twopiecewise_2050_quad', assert(3) keep(3) nogenerate
noi logout, save(output/TableS8) word dec(1) replace: tabstat yield2020 cc2050_lin_245 cc2050_quad_245 trends2050_lin_245 trends2050_quad_245 cc2050_lin_585 cc2050_quad_585 trends2050_lin_585 trends2050_quad_585, stat(p10 p25 p50 p75 p90) format(%10.4f) c(s)
save data/figure/twopiecewise_2050, replace

}

*** EOF
